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Abstract 

In CLEAN (Cryogenic Low Energy Astrophysics with Noble gases), a proposed 
neutrino and dark matter detector, background discrimination is possible if one 
can determine the location of an ionizing radiation event with high accuracy. Here, 
we develop spatial methods for event reconstruction, and study their performance 
in computer simulation experiments. We simulate ionizing radiation events that 
produce multiple scintillation photons within a spherical detection volume filled 
with liquid neon. We estimate the radial location of a particular ionizing radiation 
event based on the observed count data corresponding to that event. The count 
data are collected by detectors mounted at the spherical boundary of the detection 
volume. We neglect absorption, but account for Rayleigh scattering. To account for 
wavelength-shifting of the scintillation light, we assume that photons are absorbed 
and re-emitted at the detectors. In our study, the detectors incompletely cover the 
surface area of the sphere. In the first method, we estimate the radial location of 
the event by maximizing the approximate Poisson likelihood of the observed count 
data. To correct for scattering and wavelength-shifting, we adjust this estimate 
using a polynomial calibration model. In the second method, we predict the radial 
location of the event as a polynomial function of the magnitude of the centroid of 
the observed count data. The polynomial calibration models are constructed from 
calibration (training) data. In general, the Maximum Likelihood method estimate is 
more accurate than that of the centroid method estimate. We estimate the expected 
number of photons emitted by the event by a Maximum Likelihood method and a 
simple method based on the ratio of the number of detected photons and a detection 
probability factor. 
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1 Introduction 

We estimate the location of an ionizing radiation event that produces multiple 
photons within a spherical detection volume based on count data recorded by 
detectors mounted on the boundary of this detection volume. The detectors 
cover approximately 75 % of the total area of the detection volume boundary. 
We assume that all detected photons produced by a particular event can be 
distinguished according to their arrival times from photon counts produced by 
other events. Beyond this, we do not require additional temporal information. 
We assume that absorption is negligible, and that photons undergo Rayleigh 
scattering. We also estimate the expected number of emitted photons for the 
event. 

Our "event reconstruction" work is motivated by a proposed experiment called 
Cryogenic Low Energy Astrophysics with Noble gases (CLEAN). In CLEAN 
[1,2], events would be detected based on scintillation light produced by neutrino- 
electron scattering and WIMP-nuclear elastic scattering in a large cryostat 
filled with liquid neon. In CLEAN, the expected number of scintillation pho- 
tons produced by an event would be proportional to the amount of energy de- 
posited by the event. The proportionality factor would be determined from a 
calibration experiment. Such events of interest would occur uniformly through- 
out the cryostat. Here, we consider a cryostat with a spherical geometry. Be- 
cause neon has lower binding energy to surfaces than most radioactive impu- 
rities, there should not be internal backgrounds in the neon provided that the 
neon is purified using cold traps. Background gamma ray events would also 
produce scintillation light, and tend to occur near the walls. Near the center 
of a spherical detection volume, the penetration probability of background 
gamma rays is very low. Thus, if the radial position of an event can be deter- 
mined accurately, background gamma ray events can be discriminated with 
high confidence from events of interest. 

In this current paper, we focus solely on the event reconstruction problem. In a 
forthcoming paper, we will present a detailed simulation of CLEAN, including 
background gamma ray propagation. In this later paper, we will quantify the 
performance of our event reconstruction method for determining the location 
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of an event in the context of the physics goals of a CLEAN experiment. We 
will also discuss experimental design, calibration, and construction of CLEAN. 

Temporal methods, not spatial methods, are commonly used for event recon- 
struction in the current neutrino experiments KamLAND [3,4] and Borexino 
[5]. Both of these experiments use organic scintillators. Here, we present a 
spatial method for event reconstruction and carefully study its performance 
in computer simulation experiments. We remark that the XMASS [6] research 
team is also studying spatial event reconstruction methods for their proposed 
neutrino experiment. Liquid xenon would be used as a scintillator in XMASS. 
In the field of tomography, spatial methods are commonly used to recon- 
struct spatially varying intensity images [7-11]. For each pixel or voxel, the 
corresponding intensity parameter is the expected number of photons emitted 
during the experiment. Instead of reconstructing a spatially varying photon 
intensity image, we seek to determine the intensity and location of a single 
source of photons. Thus, our problem is a special case of a more general imag- 
ing problem. In this work, we model the observed count data as a realization 
of a spatial Poisson process. This approach is frequently taken to solve a va- 
riety of imaging problems in nuclear medicine [8-11]. We estimate the event 
location by maximizing the approximate likelihood function of the observed 
detector count data assuming a scatter-free transition matrix. Because the 
transition matrix does not account for scattering, the estimate of the radial 
location of the event can be strongly biased. However, we significantly reduce 
this bias, that is, systematic error, by correcting the initial estimate based 
on a calibration model determined from calibration (training) data. This is 
possible because the bias of the estimate introduced by misspecification of the 
transition matrix has a very predictable structure. 

We do not mean to imply that bias results only because of misspecification of 
the transition matrix. For the case where we use the exact transition matrix, 
we expect some bias because nonlinear estimation methods, including the 
Maximum Likelihood (ML) method, generally produce biased estimates. We 
expect this bias to be more significant for low count situations than for high 
count situations. For a general discussion of this point, see [12]. For a specific 
illustration of this point, see [13]. 

We stress that the methods we develop here are appropriate for the case 
where absorption of scintillation photons within the neon is negligible. For 
the case where absorption is significant, our methods are not directly appli- 
cable without further modification. If absorption effects are significant in the 
actual CLEAN experiment, we plan to directly estimate a transition matrix 
that will account for scattering, absorption, and other geometric effects. (This 
direct transition matrix modeling approach has been suggested in [6] for the 
planned XMASS experiment.) Even if we have exact knowledge of the true 
transition matrix, further calibration experiments may be necessary to quan- 
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tify the performance of CLEAN. The need for further calibration studies may 
be most pressing for low count situations. 

In Section 2, we present the scattering model. In Section 3, we present a 
Poisson likelihood model for the data based on a scatter-free transition matrix. 
In Section 4, we study the statistical properties of both the uncorrected and 
corrected radial estimates for a variety of cases. For the same data, we compare 
a prediction model based on the centroid of the observed counts with the ML 
estimate. We also estimate the expected number of emitted photons at the 
event location using the ML method and a simpler method. In the simpler 
method, the estimate is proportional to the number of detected photons. 



2 Simulation Model 

At a particular location, multiple scintillation photons are produced by an 
event of interest. The probability density function (pdf) for the initial veloc- 
ity direction of each emitted photon is uniform on the surface of the unit 
sphere. Following standard practice [14,15] we simulate the distance traveled 
before first scattering (as well as the distance between subsequent scatter- 
ings) by sampling from an exponential distribution. The expected value of 
this realization is the scattering length X s . At the point corresponding to the 
first scattering, we rotate the velocity direction about its original direction 
according to a model for the differential cross section of the Rayleigh scat- 
tering process. We neglect angular variation of the atomic form factor. Thus, 
the inner product of the new velocity direction and the old velocity direction, 
cos(0), has a pdf proportional to ( 1 + cos 2 (#) ). We select a position on 
the cone defined by cos(#) and the original velocity direction by sampling an 
azimuthal angle that is uniformly distributed between and 2n. 

We compute the point where a photon crosses the spherical boundary of the 
detection volume. If this point is not within the area of a detector, the photon is 
not detected. For the case where the crossing point is within the area covered 
by a detector element, we consider two detection models. In the "no-shift" 
model, the photon is detected by the detector element with probability p e . 
We assume that p e is independent of the photon trajectory with respect to 
the detector. In the "shift" model, the detector element absorbs the photon 
and converts it to lower energy. This wavelength-shifted photon is randomly 
re-emitted. We study the "wavelength-shifting" detection model because the 
extreme ultraviolet scintillation light that would be produced by liquid neon 
in CLEAN is not directly detectable by conventional detectors. However, after 
being shifted to a longer wavelength, this light is detectable by conventional 
detectors. The corresponding pdf of the velocity direction of the shifted photon 
is uniform on the surface of the unit sphere. If the shifted photon's velocity 
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points radially outward, the re-emitted photon is detected with probability p e 
at the detector element where shifting occurs. Otherwise, the photon travels 
inward, without scattering, until it crosses the spherical detection boundary. 
If this crossing point falls within the area covered by a detector element, it 
is detected with probability p e . Otherwise, the photon is absorbed by the 
spherical wall and lost. Our assumption that the wavelength shifted photon 
does not scatter is based on the fact that the shifted photon has much longer 
wavelength (about 400 nm) than the original scintillation photon (about 80 
nm), and Rayleigh scattering is much stronger at shorter wavelengths. 



3 Estimation Method 



3. 1 Likelihood Model 



Given that a photon is emitted at location x s , the probability that the /cth 
detector will detect this photon is denoted by P(/c|x s ,p e ) where p e is the 
efficiency of the detector. By varying k and x s , we call the set P(/c|x s ,p e ) the 
probability transition matrix. For the "no-shift" detection model where there 
is no scattering, that is, \ s = oo, the transition matrix is 



P{k\^ Pe ) = ^ I X Xs (1) 

47T J X — X s r X 



where we integrate over the area of the /cth detector. The expected number of 
detected photons at the /cth detector is 



< n k > = A fc (x s ) = XP(k\x a ,p e ), (2) 

where the intensity parameter A is the expected number of photons emitted 
during the experiment. We approximate the integral in Eq. 1 as 



P(/c|x s ,p e ) « e c [i cos(d fc -x s ,d fe ), (3) 

N det \d k - x s | 2 

where d^ is the point where the /cth detector is tangent to the sphere, p c is the 
fractional area of the surface of the sphere covered by all detectors, and N^et is 
the number of detectors. We simulate data for N det = 2072 detectors tangent 
to the spherical boundary of the detection volume (Figure 1). We determine 
the tangent points by an optimal packing scheme [16]. 
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The actual number of emitted photons at the event location is a realization 
of a Poisson process with intensity parameter A. The expected value of the 
realization equals the intensity parameter. We assume that the number of 
counts at the kth. detector is a realization of a Poisson process with intensity 
parameter X k . Thus, the Poisson log-likelihood function of the observed data 
is 

N det 

logL = -^fc(x s ) + n fc log(A fc (x s )) -log(n fc !), (4) 

k 

where n k is the number of detections at the kth detector. We estimate the 
model parameter vector 6 = (x, y, z, A) by maximizing Eq. 4. Since, we use 
an inexact transition matrix, strictly speaking, we obtain an approximate ML 
estimate of the model parameters. 

For the case where the transition matrix is exact, the (i,j)th component of 
the asymptotic covariance matrix of the ML model parameter is 



COVipJj) = (/-%, (5) 
where the Information matrix is 



I = ^ 1 dA fc (x s ) dA fc (x s ) (g) 
k=l Afc d6i d6j 

For the Poisson model, VAR(n k ) = A fc (x s ). The asymptotic variance of the kth 
parameter estimate is the kth diagonal element of the asymptotic covariance 
matrix. The asymptotic standard error (s.e.) of the k parameter estimate is 
the square root of this asymptotic variance. 

Our estimate of the radial location of the event is 



f = ^Jx 2 + y 2 + z 2 . (7) 

For the special case where x s = (0, 0, z) and there is no scattering, we 
approximate the standard error of f as 



<7 f «<7 5 = (/- 1 ) 3 3- (8) 

For more details about asymptotic statistical theory, and a discussion of a 
closely related ML estimation problem involving Poisson data, see [17] and 
[13]. 



6 



3.1.1 Efficient Simulation 



If the detector efficiency p e is much less than 1, many simulated photons 
are not detected even if they hit a detector element. Below, we describe an 
efficient simulation model which simulates the relevant fraction of emitted 
photons that are detected (provided that they hit a detector in the "no-shift" 
detection model or in the "shift-model" after wavelength shifting). On average, 
the relevant fraction that we simulate is p e . 

In our simulation model, we assume that the efficiency of the detector p e 
is independent of the photon trajectory with respect to the detector. Thus, 
AP(A;|x s ,p e ) = Ap e P(/c|x s , 1) (Eq. 2). Because of this relationship, we need not 
simulate trajectories for all simulated photons that are emitted at x s . In our 
efficient simulation model, the number of photons simulated is a realization of a 
Poisson process with expected value Xp e rather than A. However, in the efficient 
simulation model, the detector efficiency is 1 rather than p e . That is, for data 
simulated with the efficient approach, we replace P(k\x. s ,p e ) with P(/c|x s , 1). 
In both simulation models, the simulated number of detected photons at the 
kth detector is a realization of a Poisson process with intensity parameter A&. 
Because of this substitution, we estimate a modified intensity Xp e rather than 
the acutal intensity A. Throughout this work, we use this efficient approach. 

3.1.2 Optimization Details 

We seek the values of the model parameters that yield the global maximum 
(assuming that it is unique) of the approximate log-likelihood function, logL 
(Eq. 4), using an iterative algorithm. Since optimization codes search for the 
global minimum of a cost function, we actually minimize — logL. In each 
iteration, there are 3 steps. 

1. Select initial values of parameter estimates. 

2. Estimate model parameters by minimizing — logL subject to the constraint 
that the event location is within the spherical detection volume. 

3. Search for the global minimum of — logL using the estimates from step 2 
as initial values. 

In step 1, for the first iteration, the initial estimate of the location of the 
event is the centroid of the observed count data. The initial estimate of the 
modified intensity parameter \p e is the number of observed counts divided 
by the coverage factor (p c )- For subsequent iterations, the initial location is 
the sum of the centroid and a random perturbation vector. Each component 
is a normal (Gaussian) random variable with expected value equal to and 
standard deviation equal to 0.05 R. We require that the location vector has 
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magnitude less than or equal to 0.95 R. We perturb the initial estimate of Xp e 
by multiplying n/p c by a factor 1 + e, where e is a normal random variable 
with expected value and standard deviation 0.02. 

In step 2, we minimize the sum of — logL and a penalty function, using an 
algorithm based on a quasi-Newton method [18]. The penalty function forces 
the estimate to fall within the detection volume. The penalty function is 
when f < R, and 10 000 (f/R) 2 when f > R. 

In step 3, we minimize the cost function, — logL, using a modified Newton 
method where we supply both the first and second derivatives of the cost 
function with respect to the model parameters [19]. 

In general, the numerical stability of the estimate depends on the number 
of counts, the scattering length, and the geometry of the detectors. For low 
count situations, stability is usually achieved by 30 iterations. For instance, 
for the case of wavelength-shifting and X s /R = 0.1, we simulate 500 events 
of interest for each of three cases. In the three cases, the number of observed 
counts n is 10, 25 and 50. For n = 10, the difference between the uncorrected 
estimate of f/R at 30 and 90 iterations was greater than 0.01 only three 
times. For n =25 and n =50, the difference between the estimate at 30 and 
90 iterations was less than 0.01 for all 500 realizations. The number of cases 
for which the difference between the estimate at the first and 90th iterations 
was greater than 0.01 is respectively 14, 1, and 4 for n = 10, n = 25, and n = 
50. For all cases where the estimate at the first and 90th iterations varied by 
more than 0.01, the "instability" occurs at f/R > 0.65. Hence, the outcome 
of a classification rule to determine if an event is either "background" or an 
"event of interest" is less sensitive to the number of iterations than is the 
point estimate of radial location. For all cases studied here, stable estimates 
are obtained by 30 iterations, provided that the number of detected photons 
is greater than 50. 



3.2 Calibration Model 



Our calibration model is 



_£ = ai _ + a 2 (-) 2 + a,(-f, (9) 

where f c is the corrected estimate. We estimate the calibration model param- 
eters, ai,a 2 ,a 3 , from calibration data by minimizing 
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S> - «l| + « 2 (|) 2 + «3(|) 3 ) 2 , (10) 

where and f\ are the true and estimated values at the ith point in the cali- 
bration data set. The radial locations of the points in the calibration data are 
uniformly distributed (Table 1). We quantify the accuracy of the calibration 

model by computing RMSE = \j^z^YJi=\^R T r) 2 -> where m is the total 

number of points in the calibration data set and k is the number of calibration 
model parameters. For the cases studied, m ~ 2500. In Table 1 and elsewhere, 
we represent the standard error, that is, the estimated standard deviation, of 
an estimate in parentheses. For instance, 0.894(16) means the estimated value 
is 0.894, and the associated standard error is 0.016. 

In this work, we focus on estimation of the radial location of an event. If 
one wishes to estimate the Cartesian coordinates of the event, we suggest 
that the estimates of x,y and z are corrected in the same way as f. That is, 
x c = Vc = ^y, and z c = ^fz. In Figure 2, we illustrate how this approach 
works for an example. 



4 Simulation Results 

4-1 Position Estimation 

In our first study, we simulate data at selected points on the z-axis of our 
detector coordinate system for the "no-shift" detection model. For each point, 
we simulate 100 data sets. The modified intensity parameter is \p e = 200. 
For the "no-shift" case, the bias of the uncorrected estimate increases as the 
scattering length decreases (Figure 3). This bias is expected because, in gen- 
eral, scattering increases the probability that a photon emitted in the upper 
hemisphere z > will cross the spherical boundary in the upper hemisphere. 
For instance, for the case where photons are emitted at x = (0,0,0.9/?), the 
probability that the crossing occurs at z > 0.9R is 0.541(5), 0.708(5), 0.801(4) 
for X s /R = 1, 0.1, 0.01. This phenomenon appears to be related to the behavior 
of particles undergoing Brownian motion near absorbing walls [20]. 

The bias of the corrected estimate is significantly lower than the bias of the un- 
corrected estimate. In general, the RMSE of the corrected estimate decreases 
as the scattering length decreases. Moreover, the standard error of the cor- 
rected estimate is often below the asymptotic standard error computed from 
asymptotic theory for the scatter-free case. This is not inconsistent with statis- 
tical theory. Moreover, our estimation method exploits additional information 
provided in the calibration data that is not accounted for in the asymptotic 
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theory calculation. We stress that our results correspond to the case where ab- 
sorption is negligible. For the more general case where absorption effects are 
significant, the accuracy of a corresponding prediction model may decrease as 
the scattering length decreases. 

For the case where there is no scattering, that is, X s = oo, "shifting" introduces 
a negative bias (Figure 4). For the "shift" model, as the scattering length 
decreases, the bias increases in general. 

In CLEAN, events of interest would occur uniformly throughout the detection 
volume. Hence, the pdf for the radial location of an event of interest is pro- 
portional to r 2 . We simulate events of interest that produce a fixed number 
of counts n. In Figure 5 we display scatterplots of the true value of r and the 
uncorrected estimated value of r for three cases. In CASE A, we consider the 
"shift" detection model where X s /R = oo. In CASE B, we consider the "no- 
shift" detection model where X s /R = 0.1. In CASE C, we consider the "shift" 
detection model where X s /R = 0.1. In Figure 6, we display scatterplots of the 
true value of r and the corrected estimated value of r for the same three cases. 
In each scatterplot, we show the 0.1 quantile of the empirical distribution of 
the estimate as a dashed line and the level r v j R = p 1 ^ 3 where p = 0.1 as a 
solid line. If our estimate f were uniformly distributed in the spherical detec- 
tion region, the fraction of estimates less than r p = Rp 1 / 3 would be p. Thus, 
ideally, for each value of n, one tenth of the estimates should fall below r .i. 

Assuming that events that produce a fixed number of observed photons occur 
uniformly throughout the detection volume, the expected value of the frac- 
tion of corrected estimates of r below r p divided by p is the efficiency of the 
detector. Ideally, the detector efficiency should be 1 for all n. However, our 
calibration model (Eq. 9) does not guarantee a detector efficiency equal to 1. 
Since a nonuniform detector efficiency as a function of n will distort the count 
spectrum of the detected events, one should quantify the detector efficiency. 
As an illustration, we estimate detector efficiency for the "shift" detection 
model where X s /R = 2/9 (Figure 7). 

For experiments like CLEAN, we wish to determine whether an event occurs 
within a particular inner spherical region or fiducial volume defined by r < 
Rp 1 / 3 , where < p < 1. In principle, given enough calibration data, we can 
determine the radial boundary of a 100 x p percent fiducial volume so that the 
detector efficiency is arbitrarily close to 1 by setting r p equal to the pth quantile 
of the distribution of estimated values of r. This classification problem can be 
viewed as a statistical hypothesis test where the null hypothesis is that the 
actual radial location of a particular event is a realization from a distribution 
that has a pdf proportional to r 2 . In this way of looking at the problem, 
the estimated value of r (corrected or uncorrected) can be interpreted as the 
test statistic. We accept the null hypothesis if the test statistic falls below a 
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certain critical level. The pth quantile of the distribution of this test statistic 
(simulated under the assumption that the null hypothesis is true) is the critical 
level for a test with size 1 — p [21]. For a test with size 1 — p, when the null 
hypothesis is true, we reject it with probability 1 — p. When one sets the 
critical level of a test with size 1 — p to the pth. quantile of a test statistic, 
one has calibrated a confidence point [22 page 263]. We wish to reject the null 
hypothesis with high probability when an event is due to background. That is, 
we wish the power of the test to be high. When events are due to background, 
a high rejection rate will yield a low background. 



4-1.1 Centroid Method 

For comparison, we predict the radial location of an event based on the mag- 
nitude of the centroid of the count data. The centroid C is defined to be 



C = 5> fc d fc , (11) 

k 



where 



w k = Ndet (12) 

and rik is the number of counts at the kth detector. Our Centroid method 
estimate of the radial location of the event is 



f 'centroid a |C| , n (\C\\2 

— R — ~~ ~R" ^~R* ' ^ ' 

In general, the variability of the Centroid method estimate (Figure 8, Table 
2) is higher than that of the ML method (Figure 6, Table 1). The difference 
in performance is most dramatic for low count (n < 50) cases. 



4-2 Intensity Estimation 



In the efficient simulation model, we simulate a Poisson random variable N 
with expected value equal to modified intensity \p e . To derive a simple esti- 
mate of the modified intensity, we assume that the number of detected counts 
n, given N, is a binomial random variable with expected value Np and vari- 
ance Np(l — p). For the "no-shift" detection model p = p c . For the "shift" 
detection model, p = p c {l — 0.5(1 — p c )) . Thus, our "naive" estimate of the 
modified intensity is 
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Xp, 



n 



(14) 



e naive 



P 



Our second estimate of the modified intensity is 



^PeML — — #4, 



(15) 



where #4 is the ML method estimate of Xp e determined by maximizing the log- 
likelihood of the observed data (Eq. 4). We adjust 84 because our transition 



matrix does not account for the effects of wavelength-shifting by the detectors. 
For the cases studied, for r/R < 0.9, the two methods yield similar results. 
However, for r « R, the variability of ML method estimate is generally higher 
than that of the naive method estimate. To illustrate this point, we compare 
the two methods for the "shift" detector model where \ S /R = 0.1 (Figure 9). 

In an additional study, we simulate realizations of iV that are uniformly dis- 
tributed between 50 and 1200. The Cartesian coordinates of each event loca- 
tion are randomly generated so that the radial location of a random event is 
uniformly distributed between and R. We compute the test statistic 



for each event. Assuming that our binomial assumption is correct, the pdf 
of z should be well approximated by a unit normal (Gaussian) distribution 
centered on 0. That is, the expected value and variance of z should be ap- 
proximately and 1. Clearly, for events near the wall, the observed variability 
of z is inconsistent with the binomial assumption (Figure 10). We attribute 
this high variability to the incomplete coverage of the walls by the detectors. 
For instance, in some cases events can occur near the wall in gaps between 
detectors. For such cases, almost all of the emitted photons traveling radially 
outward could be lost to the walls. Assuming that the other half of the pho- 
tons that travel radially inward are detected with probability p, the overall 
detection probability would be approximately p/2 rather than p. We attribute 
the apparent outliers for CASE A in Figures 5,6, and 7 to a similar wall effect. 

For the "no-shift" detection model, the probability of detecting a photon that 
is emitted at the center of the spherical detection volume is p c = 0.75295(13). 
Based on p c , we predict the probability of detecting a photon emitted from 
the center of the spherical detection volume to be p = 0.65994(16) using 
the formula p = p c (l — 0.5(1 — p c )) presented in Section 4.2. For simulated 
events where the initial radial position of the photon was uniformly distributed 
between and 0.9 R, we compute the ratio of detected to emitted photons 



n — Np 



(16) 



z = 
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for for scattering lengths \ S /R = 0.01, 0.1, 1, oo. For the "no-shift" detection 
model, this ratio is respectively 0.7536(3), 0.7526(3), 0.7529(3), and 0.7528(3). 
For the "shift" detection model, this ratio is respectively 0.6607(3), 0.6605(3), 
0.6606(3) and 0.6604(3). 



5 Discussion 

Since the variability of the Centroid method estimate of radial location is 
higher than the variability of the ML method estimate, we expect the ML 
method to be a more powerful method for background discrimination applica- 
tions. The relative performance of the two methods should be most dramatic 
for low count situations. Quantifying this relative performance is a subject for 
further study. 

For an actual energy spectrum of interest, due to wall effects, the actual dis- 
tribution of events of interest that produce n observed events might deviate 
slightly from a uniform distribution even if the actual distribution of events 
that yield N has a uniform distribution. If this deviation from uniformity is 
significant, we should account for it when constructing critical levels of statis- 
tical tests to determine whether an event is due to a background event. 

It is reasonable to assume that the number of emitted scintillation photons is 
a realization of a Poisson process with a rate parameter proportional to the 
energy deposited by the event. Imagine that some physical process of interest 
produces events with energy deposit spectrum f(E). The spectrum f(E) is 
a pdf in energy space. Due to the Poisson assumption, the number spectrum 
g(N) generated by this process will be proportional to the convolution of f(E) 
with a transition matrix p(N\E). The spectrum for observed counts q(n) is 



where h(n\N, x s ) is a transition matrix that accounts for incomplete detector 
coverage and position effects like the wall effects noted in this work. Estimation 
of the true energy spectrum of interest f(E) from the empirical count distri- 
bution observed in an experiment is a challenging statistical inverse problem 
worthy of further study. However, it is possible to develop statistical tests to 
determine whether a theory explains the observed data without "solving" the 
inverse problem if, given a theoretical estimate of f(E), we can predict the 
count spectrum for the observed data. 



q(n) oc 




(17) 
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6 Summary 



We estimated the radial location of an ionizing radiation event that produced 
multiple photons at a point within a spherical detection volume, based on the 
observed count data collected by detectors mounted at the spherical boundary 
of the detection volume. We neglected absorption but accounted for Rayleigh 
scattering as well as photon conversion and re-emission at the detectors. In 
one method, the predicted value was a polynomial function of the centroid 
of the observed count data. In the second method, the predicted value was 
a polynomial function of the approximate maximum likelihood estimate of 
the radial location. In both cases, the polynomial correction models were con- 
structed from calibration data. In general, the ML method estimate was more 
accurate than the Centroid method estimate. In this work, our ML estimate of 
location was determined using a transition matrix that neglects scattering and 
detector conversion of photons. We corrected our estimate of radial location 
in a post-processing step using a calibration model. For the case of negligible 
absorption, as the scattering length decreased, the accuracy of the corrected 
estimate improved. For the case where absorption is significant, we may not 
observe this sort of improvement. 

The event reconstruction methods presented in this work do not account for 
photon absorption in neon. For the case where absorption is significant, we 
would either modify our methods or adopt a different method. For the case 
where absorption is significant, a spatial ML method based using estimate 
of the true transition matrix, as suggested in [6], is a promising approach. 
However, as we mentioned earlier, even if one has a perfect estimate of the 
transition matrix accounting for scattering and other effects such as absorp- 
tion, the ML method estimate of radial position may still be biased. Additional 
calibration experiments may be necessary to quantify detector efficiency even 
if we know the true transition matrix. 

We estimated the intensity of the event using the ML method and a much 
simpler scheme. In the simpler scheme, the estimate is the number of de- 
tected photons divided by a detection probability factor. In general, the sim- 
pler method was as good as or better than the ML method. For the events near 
the wall, the variability of the ML estimate was higher than the variability of 
the simple method. 

In future work, we will study how well our event reconstruction method dis- 
criminates events of interest from background events for simulated data. In 
particular, we plan to study the relative performance of the Centroid and 
Maximum Likelihood methods. We expect the relative performance of the 
two methods to be most dramatic for low count situations. For high count 
situations, both methods could be useful. For instance, the Centroid method 
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might serve to validate results obtained by the Maximum Likelihood method. 
In an actual experiment, multiple-point background events can occur because 
background gamma rays can deposit fractional amounts of their total energy. 
Additional simulation studies indicate that we can discriminate such multiple- 
point background events from single-point events of interest. 

Our work demonstrates the feasibility of spatial methods for event reconstruc- 
tion for the case where scattering is significant but absorption is negligible. 
We developed calibration models to adjust estimates of the radial location of 
an event based on simulated calibration data. In a simulation, one can sample 
calibration data throughout the entire detection volume. In a real experiment, 
one might not have such freedom. If sampled points are not sufficiently repre- 
sentative of all points in the detection volume, the accuracy of the calibration 
model could be affected. Thus, validation of empirical calibration methods for 
an actual experiment is an important topic for further study. Other important 
topics for further study include: energy spectrum estimation, quantification of 
detector efficiency and false detection rate as a function of the number of 
detected counts, and development of statistical hypothesis tests for neutrino 
models. 
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Table 1. Calibration model parameters for Maximum Likelihood Estimate. 
RMSE is the square root of the mean square prediction error, adjusted for 
degrees of freedom, computed from all calibration data. The number of 
counts in each calibration data set varies from 50 to 1200. 
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no-shift 
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Table 2. Calibration model parameters for Centroid Estimate. RMSE is the 

square root of the mean square prediction error, adjusted for degrees of 
freedom, computed from all calibration data. The number of counts in each 
training data set varies from 50 to 1200. 

X s /R ft ft RMSE 



no-shift 

00 1.483(10) 

1 1.326(6) 
0.1 1.049(5) 
0.01 0.991(4) 
shift 

00 2.094(21) 

1 1.826(15) 
0.1 1.475(8) 
0.01 1.386(7) 



0.0108(194) 0.0465 

-0.0525(115) 0.0370 

0.0167(64) 0.0291 

0.0176(53) 0.0257 

-0.113(56) 0.0675 

-0.084(35) 0.0563 

0.006(14) 0.0438 

0.018(14) 0.0415 
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Figure 1. Detector geometry for simulation study. 



20 



o 
o 



LO 

o 



(Dp 

LO 
O 



• • » 




'-0.10 -0.05 0.0 0.05 0.10 
x est/R 



o 

CD 
O 
LO 
LO 



o 

LO 

iiLO 

co-<3- 

NO 



LO 

CO 



o 

CO 



• - 1 

• • • 
• 


• 

• • • 

• 







-0.10 -0.05 0.0 0.05 0.10 
x est/R 



o 

CD 
O 

LO 
LO 

d 
o 

LO 
CO 1 * 

Op 

•-o 
d 

LO 

co 



o 
co 



»< 



-0.10 -0.05 0.0 0.05 0.10 
x est/R 



o 
d 



.To 

1° 
o 

op 



o 
d 





• 




• • 

* • 


• • • 




• 

• . *»• 




— ^* 


«». •• — 






• 


.. • • 
* • 


• 
• 





o 

CD 

d 

LO 
LO 

cc<=> 

COLO 

NLO 
0)0 

OLO 

co 

d 

o 
co 





• 

%" ••••• . 
,v • • • 


• • • 

• 


,*•••. 

• 



o 

CD 

d 

LO 
LO 

cr°' 

COLO 
>-LO 

<D ' 

bo 

OlO 

co 



o 

CO 





« «••• • 


• 


■ • . 

• 



'-0.10 -0.05 0.0 0.05 0.10 
corrected x est/R 



-0.10 -0.05 0.0 0.05 0.10 
corrected x est/R 



-0.10 -0.05 0.0 0.05 
corrected x est/R 



0.10 



Figure 2. Top row: 100 realizations of the uncorrected x,y,z and radial compo- 
nents of the Maximum Likelihood estimate of position. Bottom row: Corrected 
version of top. True position is (x,y,z) = (0,0,0.4R). The scattering length is 
X s — 0.1 R. "No-shift" detection model. The modified intensity parameter is 
\p e = 400. 
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Figure 3. Bias and standard error of corrected and uncorrections radial esti- 
mates of position for data simulated on a grid along the z-axis. The modified 
intensity \p e is 200. "No-shift" detection model. 
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Figure 4. Bias and standard error of corrected and uncorrections radial esti- 
mates of position for data simulated on a grid along the z-axis. The modified 
intensity \p e is 200. "Shift" detection model. 
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Figure 5. Uncorrected Maximum Likelihood model estimates of radial position 
for events that occur uniformly within the sphere. CASE A: As = oo and 
"shift" detection model CASE B: X s = 0.1R and "no-shift" detection model. 
CASE C: As = 0.1R and "shift" detection model. For each plot, we show the 
0.1 quantile of the prediction (dashed line) and the r p level. 
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Figure 6. Corrected Maximum Likelihood model estimates of radial position 
for events that occur uniformly within the sphere. Same data as in Figure 5. 
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Figure 7. Detector efficiency for 1 % and 10 % fiducial volumes for \ S /R = 
2/9 (± 1-sigma uncertainty intervals shown). "Shift" detection model. 



26 



CASE A 10 counts 



CASE A 50 counts 



CASE A 500 counts 




CASE B 10 counts 



CASE B 50 counts 



CASE B 500 counts 




CASEC 10 counts 



CASE C 50 counts 



CASEC 500 counts 




0.0 0.2 0.4 0.6 0.8 1.0 

r/R 



0.0 0.2 0.4 0.6 0.8 1.0 
r/R 



0.0 0.2 0.4 0.6 0.8 1.0 
r/R 



Figure 8. Predicted radial estimate based on centroid. Same data as in Figure 
5. 
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Figure 9. For the "shift" detection model and X s /R = 0.1, we compare the 
mean and standard error of two estimates of modified intensity for the case 
where event locations are on the z axis. The modified intensity is \p e = 200. 
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Figure 10. "Shift" detection model and X s /R = 0.1. We predict iV based 
on the number of detected photons n using our "naive" estimation model. 
If our "naive" model for predicting modified intensity is valid, z should be 
approximately a Gaussian random variable with expected value and standard 
deviation equal to and 1. 
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